www.gusucode.com > Matlab在化学工程中的应用 > Matlab在化学工程中的应用/实用化工计算机模拟－Matlab在化学工程中的应用/Examples/Chapter 4/ChemHeatPump.m
function ChemHeatPump
% 可逆气固化学热泵制冷系统的动态模拟--考虑热容随反应进度的变化
clear all
global R
global mt ms Ms Ns Mt Ff Cps1 Cps2 Cpf Cpt n2
global d do h1 h2 e1 e2 d1 d2 Vf Vs
global Rho
global koa kod Ea Ed Ma Md
global hsw hfw hfe Asw Afw Afe Wo
global dHe dHr IAD
global Ta Pc Tfin % Ta: 环境温度
global Toc Toh
global wgr rhos Cpgr Mg
% 输入参数
% --------
% (1)反应器结构参数
d = 0.150; % 吸附器内径，m
do = 0.006; % 气体扩散器（中心孔）的直径，m
h1 = 0.100; % 反应材料的高度，m
h2 = 0.114; % 夹套有效高度，m
e1 = 0.0045; % 壁厚，m
e2 = 0.012; % 夹套的流通宽度，m
d1 = d + 2*e1; % 吸附器外径，m (d1=0.159 m)
d2 = d+2*e1+2*e2; % 反应器外壁，m (d2=0.183 m)
Vf = pi/4*(d2^2-d1^2)*h2; % 夹套中导热油的体积，m3
Vs = pi/4*(d^2-do^2)*h1; % 反应器中盐床的体积，m3
% (2)反应器物性参数
Mt = 25; % 反应器质量，kg
Cpt = 50; % 反应器热容量，Jkg-1K-1
% (3)反应（吸附）材料的参数
mt = 1.390; % kg,IMPEX的总质量,对于SrCl2IMPEX
ms = 1037.14; % g，对于SrCl2/NH3IMPEX
Cps1 = 697; % SrCl2.NH3的热容，Jm-3K-1
Cps2 = 1480; % SrCl2.8NH3的热容，Jm-3K-1
Ms = 158.526; % 反应盐SrCl2的摩尔质量，g/mol
Mg = 17; % 氨气的摩尔质量，g/mol
n2 = 7; % 对于SrCl2/NH3
Ns = ms/Ms;
wgr = 0.30; % 可膨胀石墨的重量百分比
Cpgr = 674; % 可膨胀石墨的比热,J/kgK
rhos = mt/Vs; % 反应材料的密度，kg/m3
IAD = 1;
% (4)热力学参数
dHe = 23366; % 氨的蒸发潜热，J/mol
dHr = 41431; % 反应热，J/mol
K = 273.15;
R = 8.3145;
% (5)传热参数及传热面积
hsw = 500; % W/(Km2)
hfw = 692; % W/(Km2)
hfe = 2; % W/(Km2)
Asw = pi*d*h1;
Afw = pi*d1*h1;
Afe = pi*d2*h2;
% (6)动力学参数
koa = 0.0190;
Ea = 6921;
Ma = 2.96;
kod = 0.125;
Ed = 9000;
Md = 3.02;
% (7)载热介质(导热油)物性参数：
Rho = 870; % 载热介质（导热油）的密度，kgm-3
Cpf = 2400; % 载热介质的热容，Jkg-1K-1
% (8)操作参数
Ff = 80e-6; % 载热介质的体积流量，m3s-1
Wo = Vf*Rho; % 夹套中导热油的重量，kg
Ta = 30+K; % 环境温度，K
Th = 140 + K; % 绝对温度，K
Toc = 20 + K; % 冷油温度
Toh = 140 + K; % 热油温度
Pa = 5; % 操作压力，bar
Pd = 13;
% 吸附过程模拟
% ------------
tspan = [0 9000];
y0 = [0 Toc Toc Toc];
[t,y] = ode23s(@Equations,tspan,y0,[],1,Pa,Toc);
% 吸附过程总反应进度的动态变化图
plot(t,y(:,1))
xlabel('Time (s)')
ylabel('X_a')
title('吸附过程总反应进度的动态变化')
% 吸附过程床层温度、壁温、载热流体温度的动态变化图
figure
plot(t,y(:,2),'k--',t,y(:,3),'b-',t,y(:,4),'r-.')
xlabel('Time (s)')
ylabel('Temperature (K)')
legend('T_s','T_w','T_f')
title('吸附过程床层温度、壁温、载热流体温度的动态变化')
% 考察吸附压力的影响
Pai = [2, 3, 4, 5];
n = length(Pai);
for i = 1:n
[t,y] = ode23s(@Equations,tspan,y0,[],1,Pai(i),Toc);
ti{i} = t;
Xi{i} = y(:,1);
Tsi{i} = y(:,2);
end
% 吸附压力对总反应进度X的影响图
figure
plot(ti{1},Xi{1},'r:',ti{2},Xi{2},'b-.',ti{3},Xi{3},'k-',ti{4},Xi{4},'b--')
xlabel('Time (s)')
ylabel('X_a')
legend('Pa = 2','Pa = 3','Pa = 4','Pa = 5')
title('吸附压力对总反应进度的影响')
% 吸附压力对反应床层温度Ts的影响图
figure
plot(ti{1},Tsi{1},'r:',ti{2},Tsi{2},'b-.',ti{3},Tsi{3},'k-',ti{4},Tsi{4},'b--')
xlabel('Time (s)')
ylabel('T_s')
legend('Pa = 2','Pa = 3','Pa = 4','Pa = 5')
title('吸附压力对反应床层温度的影响')
% 解吸过程模拟
% ------------
y0 = [0 Toh Toh Toh];
[t,y] = ode23s(@Equations,tspan,y0,[],-1,Pd,Toh);
% 解附过程图形输出
figure
plot(t,y(:,1))
xlabel('Time (s)')
ylabel('X_d')
title('解吸过程总反应进度的动态变化')
figure
% plot(t,y(:,2:4))
plot(t,y(:,2),'k--',t,y(:,3),'b-',t,y(:,4),'r-.')
xlabel('Time (s)')
ylabel('Temperature (K)')
legend('T_s','T_w','T_f')
title('解吸过程床层温度、壁温、载热流体温度的动态变化')
% ------------------------------------------------------------------
function dydt = Equations(t,y,IAD,Pc,Tfin)
% 函数功能：计算某时刻的反应盐(床层)、载热流体
% 和反应器壁的温度以及反应进度的动态变化
%
% 输入参数：
% t --- 时间，s
% Cps --- 反应器内IMPEX的热容量，Jm-3K-1
% Cpf --- 载热流体（加热油）的比热容，Jkg-1K-1
% Cpt --- 反应器壁的热容量，Jkg-1K-1
% Mt --- 反应器壁的质量，kg
% Ff --- 载热流体（加热油）的体积流率，m3s-1
% Rho --- 载热介质（导热油）的密度,kgm-3
% dHr --- 反应热
% hsw --- 床层与反应器壁之间的传热参数，W/（Km2）
% Asw --- 床层的传热面积，m2
% hfw --- 反应器壁面和载热流体之间的传热参数，W/K
% hfe --- 载热流体与环境之间的散热参数，W/K
% Ta --- 环境温度，K
% Pc --- 操作压力，bar
%
% 输出变量：
% Ts --- 某时刻的反应盐(床层)温度
% Tf --- 某时刻的载热流体温度
% Tw --- 某时刻的反应器壁温度
% X --- 某时刻的总反应进度
global Ns Mt Ff Cps1 Cps2 Cpf Cpt n2
global Rho hsw hfw hfe Asw Afw Afe Wo
global dHr Ta Vf Vs
X = y(1);
Ts = y(2);
Tw = y(3);
Tf = y(4);
dXdt = Kinetics(IAD,t,X,Pc,Ts);
Cps = CPX(IAD,X);
% 对反应盐(床层) --- 吸附时IAD=1, 解吸时IAD=-1
dTsdt = ( hsw*Asw*(Tw-Ts)+IAD*n2*Ns*dHr*dXdt )/(Vs*Cps);
% 对反应器壁:
dTwdt = ( hsw*Asw*(Ts-Tw)-hfw*Afw*(Tw-Tf) )/(Mt*Cpt);
% 对载热流体:
dTfdt = ( hfw*Afw*(Tw-Tf)-hfe*Afe*(Tf-Ta) )/(Wo*Cpf) ...
- Ff*Rho*Cpf*(Tf-Tfin);
dydt = [dXdt; dTsdt; dTwdt; dTfdt];
% ------------------------------------------------------------------
function dXdt = Kinetics(IAD,t,X,Pc,Ts)
global R koa Ea Ma kod Ed Md
% koa --- 吸附系数
% kod --- 解吸系数
% Ea --- 吸附活化能
% Ed --- 解吸活化能
% Ma --- 吸附幂指数
% Md --- 解吸幂指数
if IAD == 1;
ko = koa;
E = Ea;
M = Ma;
end
if IAD == -1;
ko = kod;
E = Ed;
M = Md;
end
Ar = ko * exp(-E/R./Ts);
Peq = exp( -4983.16./Ts+16.00 );
Pm = IAD*( Pc-Peq )./Peq;
if (Pm<=0) | (X<0) % Pm <= 0表示尚未发生吸附或解吸
X = 0;
dXdt = 0;
else
arg = Ar * Pm .* (M-1) .* t + 1;
dXdt = Ar .* ( 1-X ) .^M .* Pm; % dXdt = f(Ts)
end
% ------------------------------------------------------------------
function Cp = CPX(IAD,X) % 计算给定总反应进度下的体积热容
global wgr rhos Cpgr Cps1 Cps2 Mg Ms
% wgr weight percentage of inert binder
% rhos anhydrous volumetric mass, kg/m3
% Cpgr specific heat capacity of graphite, J/kgK
% Mg molar weight of gas, kg/mole
% Ms molar weight of anhydrous salt SrCl2, kg/mol
% Cp J/m3K
if (IAD == 1)
rate = X;
else
rate = 1 - X;
end
Ms1 = rhos*(1-wgr)*(1+1*Mg/Ms)*(1-rate); % SrCl.NH3
Ms2 = rhos*(1-wgr)*(1+8*Mg/Ms)*rate; % SrCl.8NH3
Cp = rhos*wgr*Cpgr + Cps1*Ms1 + Cps2*Ms2;